/f 9 / 


Center for Turbulence Research 
Annual Research Briefs 1993 



&> 


313 



Direct simulation of isothermal- 
wall supersonic channel flow 


By Gary N. Coleman 


1* Motivation, objectives, and approach 

The motivation for this work is the fact that in turbulent flows where compress- 
ibility effects are important, they are often poorly understood. A few examples 
of such flows are those associated with astrophysical phenomena and those found 
in combustion chambers, supersonic diffusers and nozzles, and over high-speed air- 
foils. For this project, we are primarily interested in compressibility effects near 
solid surfaces. Our main objective is an improved understanding of the fundamen- 
tals of compressible wall-bounded turbulence, which can in turn be used to cast 
light upon modeling concepts such as the Morkovin hypothesis and the Van Driest 
transformation (Bradshaw 1977). 

To this end, we have performed a direct numerical simulation (DNS) study of 
supersonic turbulent flow in a plane channel with constant-temperature walls. All 
of the relevant spatial and temporal scales are resolved so that no subgrid scale 
or turbulence model is necessary. The channel geometry was chosen so that finite 
Mach number effects can be isolated by comparing the present results to well- 
established incompressible channel data (Kim, Moin & Moser 1987). Here the 
fluid is assumed to be an ideal gas with constant specific heats, constant Prandtl 
number, and power-law temperature- dependent viscosity. Isothermal- wall boundary 
conditions are imposed so that a statistically stationary state may be obtained. 
The flow is driven by a uniform (in space) body force (rather than a mean pressure 
gradient) to preserve streamwise homogeneity, with the body force defined so that 
the total mass flux is constant. 

The variables are nondimensionalized by the wall temperature, the channel half- 
width, the bulk-averaged (“mixed- mean”) density, and the bulk velocity, such that 
2 J—i pdy = 1 and | pudy = 1, where the channel walls are at y = ±1. All vari- 
ables are henceforth assumed to be dimensionless, with p representing the density, 
u = u i the streamwise velocity, (x,y, z) = (xuX 2 , x$) respectively the streamwise, 
wall-normal, and span wise coordinates, and an overbar defines an average over time 
and streamwise and spanwise directions. The nondimensional governing equations 
axe: 
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where 
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The pressure is normalized by the bulk density and bulk velocity, so the ideal gas 
law is p = pT/~iM 2 . The body force term $, is nonzero only for * = 1- The purpose 
of S, the source/sink term in (3), is explained below. Equations (1) - (3) are to be 
solved numerically subject to the isothermal, no-slip boundary conditions, 


T = 1 and u = 0 at y = ±1- 


(4) 


We therefore have as relevant nondimensional parameters (t) a Mach number, 
M, based on the bulk velocity and wall sound speed; (it) a Reynolds number, Re, 
based on the bulk density, bulk velocity, channel halfwidth, and wall viscosity; (m) 
the Prandtl number, Pr; (*'«) the ratio of specific heats, 7; and (v) the viscosity 
exponent, n, where the dynamic viscosity p oc T n . These 5 parameters are used 
to define the various DNS runs. But besides choosing appropriate values for the 
“physical” parameters, we will also artificially introduce another - to allow us to 
differentiate between mean and fluctuation compressibility effects. 

The Mach number appears in the energy equation (3) in the term that represents 
the irreversible loss of kinetic energy into heat. Following Buell (1991), we interpret 
(in our simulations) the actual Mach number M in (2) and the “dissipation Mach 
n um ber” M d in (3) as separate parameters. By setting Md to values different from 
M in the DNS, we produce an effective heat source/sink S in (3) which is given by 


S = (Mj- M 2 ) 


7(7 ~ 1) r »j 
Re p dxj 


(5) 


Consequently, we can consider cases with different mean temperature profiles (and 
thus different mean property variations) at the same M. Results from these un- 
physical” M ^ Md DNS runs can, therefore, be used to determine the relative im- 
portance of turbulent-fluctuation and variable-property influences at a given Mach 

number. . _ 

Three DNS cases will be discussed, with the Mach number ranging from M — l.o 
to 3 All the runs share the same Prandtl number, specific heat ratio, and viscosity 
exponent {Pr = 0.7, 7 = 1.4 and n = 0.7), while the Reynolds number (for reasons 
given below) is either 3000 or 4880. A summary of the parameters is listed in 
table 1. Cases denoted by a single letter (A and B) in table 1 represent “physical 
simulations for which M d = M. For the M d ± M run, Case AX, M = 1.5 and 
M d = 0 Since the temperature fields in both the physical and unphysical runs 
depend almost exclusively on M d (Coleman et al. 1993), this parameter combination 
will produce the behavior of the “extra” source/sink S, eq. (5), that is necessary 
to isolate mean and fluctuation effects. With M d = 0, S is such that the mean 
temperature and density are constant across the channel, as we shall see below. 

The DNS results were generated using the code developed by Buell to study 
compressible Couette flow. During the computations, the body force <f>, is adjusted 
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Table 1. DNS physical parameters. 


Case 

M 

Mi 

Re 

Pr 

7 

M 

A 

1.5 

1.5 

3000 

0.7 

1.4 

y0.7 

B 

3 

3 

4880 

0.7 

1.4 

y0.7 

AX 

1.5 

0 

3000 

0.7 

1.4 

jiO.7 


Table 2. DNS numerical parameters. 


nx 

ny 

nz 

nxc 

nyc 

nzc 

L z 

L t 

no 

90 

60 

144 

119 

80 

4n 

4tt/3 


so that the total mass flux through the channel remains constant. (Once the flow 
reaches a statistically stationary state the variations of 4>, with time are small). The 
code utilizes a Fourier- Legendre spectral discretization along with a hybrid implicit- 
explicit third-order time-advance algorithm designed to maximize the range of Mach 
numbers that may be considered (Buell 1991). The numerical parameters used by 
all three runs are given in table 2, where L x and L z are the streamwise and spanwise 
domain sizes, and ( nx,ny,nz ) and ( nxc,nyc,nzc ) are respectively the number of 
expansion coefficients and collocation (quadrature) points in the streamwise, wall- 
normal, and spanwise directions. The runs were made on the CCF and NAS Cray 
YMP and C-90 computers at NASA Ames Research Center. 

2. Results 

An indication of the numerical fidelity of the DNS is provided by the streamwise 
and spanwise one-dimensional spectra from the channel centerline and near the 
walls, shown in figure 1. They are typical of those found from all three DNS runs in 
their rapid fall-off at high wavenumber, which implies that the x- and 2 -resolution 
is adequate. The “Legendre spectra” (not shown) also verify that the wall-normal 
resolution is sufficient. The high wavenumber streamwise and spanwise spectra 
of the velocity at both the channel centerline (figure l(a,b)) and near the wall 
(figure l(c,d)) are similar to those found in the incompressible channel (see figure 
3 of Kim et al. 1987). In the present simulations, we find that the density and 
temperature spectra are closely related to each other and that their magnitudes are 
much larger near the walls than they tire at the centerline. The streamwise and 
spanwise correlations in figure 2 are also roughly equivalent to the incompressible 
results (cf. figure 2 of Kim et al. 1987) except for two characteristics: the large 
spanwise coherence of the density and temperature at the centerline (figure 2b), 
and the greater streamwise coherence of the p, u, and T fields near the wall (figure 
2c). V 6 



FIGURE 1 . One-dimensional spectra for Case A: , 

w - , T. (a) k (c) Streamwise; (b) k (d) spanwise 
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Figure 2. Two-point correlations for Case A: symbols as in figure 1. (a) & (c) 
Streamwise; (b) & (d) spanwise. 



FIGURE 3. Two-point density correlations for Case A at channel centerline: 
full field; , with acoustic eigenfunctions removed. 








318 


G. N. Coleman 


We first discuss the spanwise coherence at the centerline, which is thought to be 
due to acoustic resonance. Evidence for this is provided first by the fact that 
the coherence is not present in the velocity, and most significantly by the re- 
sults in figure 3, which contrast the spanwise density correlation from an instan- 
taneous field (one which contributed to the figure 2b curve) with that^ obtained 
by eliminating the contribution from several “acoustic eigenfunctions. This is 
done by projecting the DNS fields on eigenfunctions of the linear inviscid lsen- 
tropic problem for a given base flow. The acoustic density and velocity fluc- 
tuations are respectively assumed to satisfy p a (x,t) — Z,k P\*-,y) e an 

u a/ x t) = Y Wi(k,y)e i(k x_w<) . At a given wavenumber k = (k x ,k z ), the lin- 

earized, isentropic Euler equations in Fourier space in terms of the “acoustic eigen- 
functions" q,(k,») = (p(k, y), nth, y), u(k, y), tti(k, y))J can then be written as 

£(q) = uq, where 


( k x up + k x pu — i(pv) y + k z pw \ 


£(q) = 


2 ^ 


V 


k x uu — i (du/dy)v + k x ^p 
k x uv — ia 2 (p/p) y 
k x uw + k z yp 


( 6 ) 


with a equal to the sound speed and v = 0aty — ±1. 

The projection of the full DNS field onto the acoustic subspace is performed by 
computing the inner product J+ 1 q DNS • q * m dy of the DNS field with eigenfunctions 

qj^ of the adjoint problem to (6), £*(q*) = ^q*> suc ^ that /_ j £(q) ■ c\ dy - 
q . £*(q *)dy. The adjoint operator is 


r*(q*) 


/k x up+ + k x ^u*+ ^(a 2 u *) 9 d -k z £wA 
k x uu * + k x pp* 

k x ud+ + i p(p+)y ~ K<te/dy)u* 

k x tiw* + k z ppi, ) 


(?) 


where u* = 0 at y = ±1. ........ ( , R s 

When the base flow is uniform (no y variation), the eigenfunctions from (b) 

are irrotational, and the eigenvalues u from (6) and (7) give phase speeds c x = 

Real(u;)/fc I that satisfy 

u — c x = ±a[(£7r/2fc x ) 2 + 1] 1/2 , ( 8 ) 

where t is the wall-normal wave number (equivalent to the number of times |p| 
and |u| change sign between -l<y<d-l). Here we use the Case A mean profiles 
shown in figure 4 as the base flow and numerically compute solutions to (6) and (7), 
which leads to “acoustic” (isentropic) eigenfunctions that do not satisfy the above 
phase relation and, in fact, have nonzero vorticity near the walls. Near the channel 
centerline, however, the eigenfunctions have a more typically acoustic behavior in 
that their ratio of dilation to enstrophy is very large and in that a given t mode (now 
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FIGURE 5. Local Reynolds number profiles: , Case A; , Case B. 

defined as the number of \p\ sign changes) consists of an upstream- and downstream- 
propagating pair* with positive and negative phase speeds relative to the centerline 
velocity, u c — c x . The dashed curve in figure 3 represents the density field after the 
isentropic modes q/ that are recognized as acoustic in the range £ = [0, . . . ,4] have 
been projected and removed, for k x L x /2n = [0, . . . , +4] (using conjugate symmetry 
to account for k x < 0) and k z L z j 2ir = [—4, . . . , +4]. Only the eigenfunctions with 
(a) very large dilation-to-enstrophy ratio near the centerline, (b) an equal number 
of sign changes for \p\ and for |Sj, and (c) no more than two modes at each £ - 
one with positive and one with negative relative phase speed - were chosen from 
the full inviscid isentropic function space to be included in the projection. The 
second criterion is overly conservative: it is useful when automating the selection 

^ Although some £ modes appear to have only a single downstream-propagating component. 
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FIGURE 6. Near-wall two-point correlations for Case B: symbols as in figure 1. 
(a) Streamwise; (b) spanwise. 

process, but excludes a few modes (because of low amplitude “wiggles” in |u|) that 
are thought to be acoustic rather than vortical. For this reason, and perhaps also 
because of physical differences between the uniform and variable mean cases that 
produce at certain l no non-vortical modes moving upstream with respect to the 
centerline velocity, the projection at some k did not include two modes for every £ in 
the 0 to 4 range. Nevertheless, the magnitude of the reduction in figure 3 suggests 
that there are significant acoustic disturbances within the simulation results. The 
correlation at z « 2 would presumably be reduced still further if more k 2 = 0 
modes were used in the projection. Note that because the computations assume 
that the channel walls are perfectly rigid (and use periodic boundary conditions), 
any acoustic signals present in the DNS are not necessarily expected to be identical 
to those found in a laboratory wind tunnel since in the simulations there is no 
mechanism for the acoustic energy to radiate away. 

The other difference, mentioned above, between the two-point correlations for the 
present and incompressible DNS is in the larger near-wall streamwise correlations 
found in figure 2c; this indicates that the near-wall streaks, which are characteristic 
of wall-bounded turbulent flows (Robinson 1991), are more coherent in Case A 
than in the incompressible channel results. At first glance, it might appear that the 
streak modification is a low Reynolds number effect (so that the effective streamwise 
domain size in wall units is smaller) since the Reynolds number based on mean 
centerline velocity is higher in the incompressible DNS than that found here. The 
variation of the local Reynolds number across the channel is shown by the dashed 
curve in figure 5, and the centerline value (2770) is seen to be slightly less than the 
3300 quoted for the incompressible channel data (Kim et al. 1987). But because the 
isothermal boundary conditions lead to a flow with a maximum mean temperature 
near the centerline and maximum density at the walls (figure 4), the mean kinematic 
viscosity ji/p (where Ji = T n ) is maximum at y = 0. Therefore, the local Reynolds 
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FIGURE 7. Turbulent Mach number, M t = M\J uJu'/VT: , Case A; 

Case B. 



FIGURE 8. Root mean square density fluctuations: , Case A; , Case B. 


number in the present DNS is apt to be larger near the walls than for when v = 
constant. This suggests that the enhanced near-wall coherence in figure 2c is solely 
a compressibility effect, a fact that is reinforced by the Case B results. In order 
to obtain a local Reynolds number profile RepuS/Ji (where 6 = 1) at M = 3 
(Case B) that remains comparable to that for Case A, the bulk Reynolds number 
was increased from 3000 to 4880. As the Case B profile (dotted curve) in figure 5 
shows, the Reynolds number is in fact at any y slightly larger than that for Case 
A (dashed), which implies that the further increase shown in figure 6a (over that 
seen in figure 2c) of the near-wall streamwise correlation for Case B is not a viscous 
effect. 

It therefore appears that the extra coherence is due to compressibility, although 
its precise source is at this point uncertain. One possibility is near-wall viscosity 
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FIGURE 10. Pressure-dilatation correlation: , Case A; , Case B. 


fluctuations (Tritton 1961; Bradshaw & Ferriss 1971); another is small-scale acous- 
tic fluctuations that are “channeled” along the cold, low-speed streaks, which act as 
an acoustic “wave-guide.” Fairly large turbulent Mach numbers and r.m.s. density 
fluctuations are found in both flows, especially near the walls (figures 7 & 8), which 
might be evidence of significant dilational effects. However, the dilational field as- 
sociated with the near-wall fluctuations is not so important els to directly increase 
the turbulent kinetic energy dissipation rate to any great degree. This can be seen 
from figure 9, which gives the ratio of the mean-square dilatation fluctuations to 
those of the mean-square vorticity: the ratio of dilatational-to-solenoidal homoge- 
neous kinetic energy dissipation (Zeman 1990; Blaisdell, Mansour & Reynolds 1993; 
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FIGURE 12. Mean velocity: (a), in wall units; (b), with Van Driest transformation; 
, M = 0 (Kim et al 1987); , Case A; , Case B; , 2.44 In + 5.2. 
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FIGURE 13. Reynolds- stress correlation coefficient: , M — 0 (Kim et 

al. 1987); , Case A; , Case B. 


Speziale & Sarkar 1991; Lele 1994). While the ratio increases by an order of magni- 
tude as M increases from 1.5 to 3, it never becomes significantly larger than 10~ 3 . 
On the other hand, the pressure-dilation correlation is found to be larger than 50% 
of the solenoidal dissipation for both Mach numbers. Figure 10 shows that near 
the walls the dilatational field creates a strong sink of turbulent kinetic energy as 
kinetic energy is transferred to the pressure fluctuations (Blaisdell et al. 1993; Lele 
1994), while toward the centerline, the pressure-dilatation acts as a smaller - but 
still important - source of kinet ic ene rgy. In the future, we hope to understand the 
link between the large negative p'uV and the observed wall-streak modification. 

With such large dilational effects present, one might surmise that this flow is 
not governed by Morkovin’s hypothesis (Favre 1992), which states that relation- 
ships between statistical properties of turbulence are unaffected by compressibility 
if the r.m.s. density fluctuations are small (of order 1/10) compared to the abso- 
lute density (Bradshaw & Ferriss 1971; Bradshaw 1977; Spina et al. 1994). But 
the density fluctuations for both Mach numbers are within the allowed range of 
0(1/10) (Figure 8), and for at lease one important statistical ratio, the mixing 
length (-uV) 1/2 /(diZ/dy), Morkovin’s hypothesis is found to work fairly well. Fig- 
ure 11 demonstrates that this quantity is reasonably independent of Mach number. 

With the invariance of the mixing length established, the so-called Van Driest 
transformation for the mean velocity immediately follows. That is, the density- 
weighted mean velocity 



(where ~p w is the mean density at the wall and the 4- superscript denotes wall units), 
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Figure 14. Mean profiles for Case AX: , u ; , ~ p \ , T. 

is expected to satisfy the incompressible log law, 

ut D = - lny + + C, (10) 

with k and C similar to their incompressible values, k « 0.4, and C « 5.2 (Bradshaw 
1977; Huang, Bradshaw & Coakley 1993; Huang & Coleman 1993). The mean 
velocity in both wall units and the Van Driest form is plotted in figure 12 (using for 
the latter a mixing length formulation for the mean temperature to write Tiy D as a 
function of T f*, the surface heat flux and the mean surface temperature (Bradshaw 
1977)). The agreement of the curves in figure 12b, especially their slopes, tends to 
reinforce the validity of the Van Driest transformation (cf. Huang Sz Coleman 1993) 
and, by extension, the Morkovin hypothesis. 

Not all statistical ratios are found to be independent of Mach number, however, 
as the Re ynol ds-stress correlation coefficients in figure 13 show. The near-wall max- 
imum of lu't/l/uj-msVrmg increases from less then 0.5 for the incompressible channel 
to over 0.6 for M = 3. (Note that for M = 0, this correlation coefficient does 
not vary appreciably with Reynolds number (Kim et al. 1987), which points to 
compressibility and not viscous effects as the source of the differences in figure 13.) 

It thus appears that the isothermal-wall channel contains some “non- Morkovin 11 
phenomena. However, it would not at this point be appropriate to firmly state that 
the results in figure 13 represent a formal contradiction to the Morkovin hypothesis 
since the hypothesis does not (regardless of the density fluctuation level) claim to 
account for the influence of spatial gradients of the mean properties (Bradshaw 
1977), which are apt to be important for this flow. Evidence of just how important 
can be found in the results from Case AX, for which Md = 0 so that the mean density 
and temperature are constant (figure 14). The near-wall streamwise correlations for 
Case AX are given in figure 15. No indication of the enhanced streak coherence 
found for Cases A and B is observed (cf. figures 2(c) & 6(a)). Therefore, wall- 
normal gradients of the mean properties are required for the streak modification to 
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FIGURE 15. Near- wall two-point streamwise correlations for Case AX: symbols as 
in figure 1. 



FIGURE 16. Pressure-dilatation correlation: , Case A; , Case B; 

Case AX. 


occur. The mean gradients are also necessary for the near-wall fluctuation effects 
(which are presumably related to the streak coherence) to be present, as is shown 
by the Case AX pressure-dilatation profile in figure 16 (solid curve). Com pared to 
the Md = 1-5 and 3 results, when the mean properties are uniform, p'u' t is much 
less important and represents a near-wall source rather than sink of kinetic energy. 
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3. Future plans 

An immediate task is to attempt to resolve the open questions regarding the extra 
streak coherence induced by the compressibility. Namely: 

A) How do the dilatation fluctuations inferred by the large pressure-dilatation cor- 
relations influence the vortical field and hence the streak structure? 

B) In what sense is the enhanced streak coherence related to the increased u'v* 
correlation coefficient? 

C) How do the variable-mean properties couple with near-wall dilatational fluctua- 
tions? 

Recommended long-term efforts include comparing the present results to those 
computed for the adiabatic- wall channel (which to develop a statistical equilibrium 
will require either one wall to be isothermal or for the flow to contain a distributed 
heat sink). It should also be of interest to compare channel results to compress- 
ible boundary layer turbulence and thus ascertain the importance of the acoustic 
disturbances that are trapped between the channel walls but free to radiate away 
in a boundary layer. New numerical schemes, possibly using a fully implicit time- 
advance algorithm, should be developed to allow study of turbulent compressible 
flows in the hypersonic regime. 
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